function [num, den, z, p] = elliptic(n, Rp, Rs, Wn, varargin)
%ELLIP Elliptic or Cauer digital and analog filter design.
%   [B,A] = ELLIP(N,Rp,Rs,Wn) designs an Nth order lowpass digital
%   elliptic filter with Rp decibels of ripple in the passband and
%   a stopband Rs decibels down.  ELLIP returns the filter
%   coefficients in length N+1 vectors B (numerator) and A (denominator).  
%   The cut-off frequency Wn must be 0.0 < Wn < 1.0, with 1.0 corresponding 
%   to half the sample rate.  Use Rp = 0.5 and Rs = 20 as starting
%   points, if you are unsure about choosing them.
%
%   If Wn is a two-element vector, Wn = [W1 W2], ELLIP returns an 
%   order 2N bandpass filter with passband  W1 < W < W2.
%   [B,A] = ELLIP(N,Rp,Rs,Wn,'high') designs a highpass filter.
%   [B,A] = ELLIP(N,Rp,Rs,Wn,'stop') is a bandstop filter if Wn = [W1 W2].
%
%   When used with three left-hand arguments, as in
%   [Z,P,K] = ELLIP(...), the zeros and poles are returned in
%   length N column vectors Z and P, and the gain in scalar K. 
%
%   When used with four left-hand arguments, as in
%   [A,B,C,D] = ELLIP(...), state-space matrices are returned. 
%
%   ELLIP(N,Rp,Rs,Wn,'s'), ELLIP(N,Rp,Rs,Wn,'high','s') and 
%   ELLIP(N,Rp,Rs,Wn,'stop','s') design analog elliptic filters.  
%   In this case, Wn can be bigger than 1.0.
%
%   See also ELLIPORD, CHEBY1, CHEBY2, BUTTER, BESSELF, FREQZ, FILTER.

%   References:
%     [1] T. W. Parks and C. S. Burrus, Digital Filter Design,
%         John Wiley & Sons, 1987, chapter 7, section 7.3.3.

[btype,analog,errStr] = iirchk(Wn,varargin{:});
error(errStr)

if n>500
	error('Filter order too large.')
end

% step 1: get analog, pre-warped frequencies
if ~analog,
	fs = 2;
	u = 2*fs*tan(pi*Wn/fs);
else
	u = Wn;
end

% step 2: convert to low-pass prototype estimate
if btype == 1	% lowpass
	Wn = u;
elseif btype == 2	% bandpass
	Bw = u(2) - u(1);
	Wn = sqrt(u(1)*u(2));	% center frequency
elseif btype == 3	% highpass
	Wn = u;
elseif btype == 4	% bandstop
	Bw = u(2) - u(1);
	Wn = sqrt(u(1)*u(2));	% center frequency
end

% step 3: Get N-th order elliptic lowpass analog prototype
[z,p,k] = ellipap(n, Rp, Rs);

% Transform to state-space
[a,b,c,d] = zp2ss(z,p,k);

% step 4: Transform to lowpass, bandpass, highpass, or bandstop of desired Wn
if btype == 1		% Lowpass
	[a,b,c,d] = lp2lp(a,b,c,d,Wn);
%
%elseif btype == 2	% Bandpass
%	[a,b,c,d] = lp2bp(a,b,c,d,Wn,Bw);

%elseif btype == 3	% Highpass
%	[a,b,c,d] = lp2hp(a,b,c,d,Wn);

%elseif btype == 4	% Bandstop
%	[a,b,c,d] = lp2bs(a,b,c,d,Wn,Bw);
end

% step 5: Use Bilinear transformation to find discrete equivalent:
if ~analog,
	[a,b,c,d] = bilinear(a,b,c,d,fs);
end

if nargout == 4
	num = a;
	den = b;
	z = c;
	p = d;
else	% nargout <= 3
% Transform to zero-pole-gain and polynomial forms:
	if nargout == 3
		[z,p,k] = ss2zp(a,b,c,d,1);
		num = z;
		den = p;
		z = k;
	else % nargout <= 2
		den = poly(a);
                [z,k] = tzero(a,b,c,d);
                num = k * poly(z);
                num = [zeros(1,length(den)-length(num))  num];
	end
end

function [z,p,k] = ellipap(n, rp, rs)
%ELLIPAP Elliptic analog lowpass filter prototype.
%   [Z,P,K] = ELLIPAP(N,Rp,Rs) returns the zeros, poles, and gain
%   of an N-th order normalized prototype elliptic analog lowpass
%   filter with Rp decibels of ripple in the passband and a
%   stopband Rs decibels down.

%   Author(s): L. Shure, 3-8-88
%   Copyright (c) 1988-98 by The MathWorks, Inc.
%   $Revision: 1.15 $  $Date: 1997/11/26 20:13:39 $

%   References:
%     [1] T. W. Parks and C. S. Burrus, Digital Filter Design,
%         John Wiley & Sons, 1987, chapter 7, section 7.3.7-8.

if n == 1,     % special case; for n == 1, reduces to chebyshev type 1
	z = [];
	p = -sqrt(1/(10^(rp/10)-1));
	k = -p;
	return
end

epsilon = sqrt(10^(0.1*rp)-1);	%rp, dB of passband ripple
k1 = epsilon/sqrt(10^(0.1*rs)-1);	%rs, dB of stopband ripple
k1p = sqrt(1-k1^2);
if k1p == 1
	error('Cannot design filter, Rp and Rs specifications are too strict.')
end
wp = 1;	% passband edge - normalized
if abs(1-k1p^2) < eps
	krat = 0;
else
	capk1 = ellipke([k1^2,k1p^2]);
	krat = n*capk1(1)/capk1(2);	% krat = K(k)/K'(k) -- need to find relevant k
end

% try to find m, elliptic parameter, so that K(m)/K'(m) = krat -- K, complete
% elliptic integral of first kind
fopt = foptions;
fopt(14) = 250;
m = fmins('kratio',.5,fopt,[],krat); 
if m<0 | m>1
    m = fmin('kratio',0,1,fopt,krat);
end

capk = ellipke(m);
ws = wp/sqrt(m);	% stopband edge (=> transition band is ws-wp in width)
m1 = 1 - m;

% find zeros; they are purely imaginary and paired in complex conjugates
j = (1-rem(n,2)):2:n-1;
[ij,jj] = size(j);
% s is Jacobi elliptic function sn(u)
[s,c,d] = ellipj(j*capk/n,m*ones(ij,jj));
is = find(abs(s) > eps);
z = 1 ./(sqrt(m)*s(is));
z = i*z(:);	% make column vector
z = [z ; conj(z)];
% order the zeros for convenience later on
z = cplxpair(z);

% poles; one purely real if n is odd - the remainder are complex conjugate pairs
% calculate v0, a 'fundamental' parameter for the poles related to inverse sc
% function . I.e. find r so sn(r)/cn(r) = 1/epsilon for the given parameter mp
r = fmins('vratio', ellipke(1-m),fopt,[],1/epsilon,k1p^2);
v0 = capk*r/(n*capk1(1));
[sv,cv,dv] = ellipj(v0,1-m);
p = -(c.*d*sv*cv + i*s*dv)./(1-(d*sv).^2);
p = p(:);   % make column vector
% check to see if there is a real pole
if rem(n,2)
	ip = find(abs(imag(p)) < eps*norm(p));
	[pm,pn] = size(p);
	pp = 1:pm;
	pp(ip) = [];
	p = [p ; conj(p(pp))];
else
	p = [p; conj(p)];
end
p = cplxpair(p);	% order poles for later use

% gain
k = real(prod(-p)/prod(-z));
if (~rem(n,2))	% n is even order so patch gain
	k = k/sqrt((1 + epsilon^2));
end

function [order,wn] = ellipord(wp,ws,rp,rs,opt)
%ELLIPORD Elliptic filter order selection.
%   [N, Wn] = ELLIPORD(Wp, Ws, Rp, Rs) returns the order N of the lowest 
%   order digital elliptic filter that loses no more than Rp dB in the
%   passband and has at least Rs dB of attenuation in the stopband.  
%   Wp and Ws are the passband and stopband edge frequencies, normalized 
%   from 0 to 1 (where 1 corresponds to pi radians). For example,
%       Lowpass:    Wp = .1,      Ws = .2
%       Highpass:   Wp = .2,      Ws = .1
%       Bandpass:   Wp = [.1 .8], Ws = [.2 .7]
%       Bandstop:   Wp = [.2 .7], Ws = [.1 .8]
%   ELLIPORD also returns Wn, the elliptic natural frequency to
%   use with ELLIP to achieve the specifications.
%
%   [N, Wn] = ELLIPORD(Wp, Ws, Rp, Rs, 's') does the computation for an 
%   analog filter, in which case Wp and Ws are in radians/second.
%
%   NOTE: If Rs is much much greater than Rp, or Wp and Ws are very close, 
%   the estimated order can be infinite due to limitations of numerical 
%   precision.
%       
%   See also ELLIP, BUTTORD, CHEB1ORD, CHEB2ORD.

%   Author(s): L. Shure, 6-9-88
%              T. Krauss, 11-18-92, updated
%   Copyright (c) 1988-98 by The MathWorks, Inc.
%   $Revision: 1.13 $  $Date: 1997/12/02 19:16:49 $

%   Reference(s):
%       [1] Rabiner and Gold, p 241.

if nargin == 4
	opt = 'z';
elseif nargin == 5
	if ~strcmp(opt,'z') & ~strcmp(opt,'s')
		error('Invalid option for final argument.');
	end
end
np1 = length(wp);
ns1 = length(ws);
if (np1 ~= ns1)
	error('The frequency vectors must both be the same length.')
end
ftype = 2*(np1 - 1);
if wp(1) < ws(1)
	ftype = ftype + 1;	% low (1) or reject (3)
else
	ftype = ftype + 2;	% high (2) or pass (4)
end

% first, prewarp frequencies from digital (unit circle) to analog (imag. axis):
if strcmp(opt,'z')	% digital
	WP=tan(pi*wp/2);
	WS=tan(pi*ws/2);
else  % don't have to if analog already
	WP=wp;
	WS=ws;
end

% next, transform to low pass prototype with passband edge of 1 and stopband
% edges determined by the following: (see Rabiner and Gold, p.258)
if ftype == 1	% low
	WA=WS/WP;
elseif ftype == 2	% high
	WA=WP/WS;
elseif ftype == 3	% stop
	fo = foptions;
	wp1 = fmin('bscost',WP(1),WS(1)-1e-12,fo,1,WP,WS,rs,rp,'ellip');
	WP(1) = wp1;
	wp2 = fmin('bscost',WS(2)+1e-12,WP(2),fo,2,WP,WS,rs,rp,'ellip');
	WP(2) = wp2;
	WA=(WS*(WP(1)-WP(2)))./(WS.^2 - WP(1)*WP(2));
elseif ftype == 4	% pass
	WA=(WS.^2 - WP(1)*WP(2))./(WS*(WP(1)-WP(2)));
end

% find the minimum order elliptic filter to meet the more demanding spec:
WA = min(abs(WA));
epsilon = sqrt(10^(0.1*rp)-1);
k1 = epsilon/sqrt(10^(0.1*rs)-1);
k = 1/WA;
capk = ellipke([k^2 1-k^2]);
capk1 = ellipke([(k1^2) 1-(k1^2)]);
order = ceil(capk(1)*capk1(2)/(capk(2)*capk1(1)));

% if both warnings are in effect, only print the first one
if (1-k1.^2) == 1
   warning('(ellipord) attenuation too strenuous, estimated order is infinite.')
elseif k^2 == 1
   warning('(ellipord) band edges too close, estimated order is infinite.')
end

% natural frequencies are simply the passband edges (WP).
% finally, transform frequencies from analog to digital if necessary:
if strcmp(opt,'z')	% digital
	wn=(2/pi)*atan(WP);  % bilinear transform
else
	wn=WP;
end

function a = kratio(m,krat)
%KRATIO Utility function for use with ELLIP.
%   KRATIO(m,krat) is a function used to calculate the zeros of an
%   elliptic filter.  It is used with FMINS to find a parameter m 
%   satisfying ellipke(m)/ellipke(1-m) = krat.

%   Copyright (c) 1988-98 by The MathWorks, Inc.
%   $Revision: 1.10 $  $Date: 1997/11/26 20:13:40 $

% to ensure we don't call ellipke(1) which is inf on non-ieee machines
% and that we only call with positive m.
m = min(1,max(m,0));
if abs(m) > eps & abs(m)+eps < 1
	k = ellipke([m,1-m]);
	r = k(1)./k(2) - krat;
elseif abs(m) <= eps	% m==0
	r = -krat;
else	% m==1 => r == inf, but can't for non-ieee machines
	r = 1e20;
end
a = abs(r);

function [a,b,c,d] = zp2ss(z,p,k)
%ZP2SS	Zero-pole to state-space conversion.
%	[A,B,C,D] = ZP2SS(Z,P,K)  calculates a state-space representation:
%		.
%		x = Ax + Bu
%		y = Cx + Du
%
%	for a system given a set of pole locations in column vector P,
%	a matrix Z with the zero locations in as many columns as there are
%	outputs, and the gains for each numerator transfer function in
%	vector K.  The A,B,C,D matrices are returned in block diagonal
%	form.
%
%   	This function handles SIMO systems if the Control Toolbox is
%   	present and SISO systems if only the Signal Processing Toolbox
%   	is installed.
%   	
%	See also SS2ZP.


[mn,nn] = size(z); [md,nd]=size(p);
% Put it in column format if its SISO and in row format.
if (length(k)==1 & md < 2 & mn < 2) & (nn > 1 | nd > 1) 
	z = z'; p = p'; 
end
[m,n] = size(z);
if n==0, n=length(k); end  % Fix to handle multi-output when z is empty
if length(k) ~= n & (~isempty(z))
	error('Z should be a column vector or K should be SIMO.');
end
if n > 1
	% If it's multi-output, we can't use the nice algorithm
	% that follows, so use the numerically unreliable method
	% of going through polynomial form, and then return.
	eval('[num,den] = zp2tf(z,p,k);') % Suppress compile-time diagnostics
	[a,b,c,d] = tf2ss(num,den);
	return
end

% Strip infinities and throw away.
p = p(finite(p));
z = z(finite(z));

% Group into complex pairs
np = length(p);
nz = length(z);
z = cplxpair(z,1e6*nz*norm(z)*eps + eps);
p = cplxpair(p,1e6*np*norm(p)*eps + eps);

% Initialize state-space matrices for running series
a=[]; b=[]; c=[]; d=1;

% If odd number of poles AND zeros, convert the pole and zero
% at the end into state-space.
%	H(s) = (s-z1)/(s-p1) = (s + num(2)) / (s + den(2))
if rem(np,2) & rem(nz,2)
	a = p(np);
	b = 1;
	c = p(np) - z(nz);
	d = 1;
	np = np - 1;
	nz = nz - 1;
end

% If odd number of poles only, convert the pole at the
% end into state-space.
%  H(s) = 1/(s-p1) = 1/(s + den(2)) 
if rem(np,2)
	a = p(np);
	b = 1;
	c = 1;
	d = 0;
	np = np - 1;
end	

% If odd number of zeros only, convert the zero at the
% end, along with a pole-pair into state-space.
%   H(s) = (s+num(2))/(s^2+den(2)s+den(3)) 
if rem(nz,2)
	num = real(poly(z(nz)));
	den = real(poly(p(np-1:np)));
	wn = sqrt(prod(abs(p(np-1:np))));
	if wn == 0, wn = 1; end
	t = diag([1 1/wn]);	% Balancing transformation
	a = t\[-den(2) -den(3); 1 0]*t;
	b = t\[1; 0];
	c = [1 num(2)]*t;
	d = 0;
	nz = nz - 1;
	np = np - 2;
end

% Now we have an even number of poles and zeros, although not 
% necessarily the same number - there may be more poles.
%   H(s) = (s^2+num(2)s+num(3))/(s^2+den(2)s+den(3))
% Loop thru rest of pairs, connecting in series to build the model.
i = 1;
while i < nz
	index = i:i+1;
	num = real(poly(z(index)));
	den = real(poly(p(index)));
	wn = sqrt(prod(abs(p(index))));
	if wn == 0, wn = 1; end
	t = diag([1 1/wn]);	% Balancing transformation
	a1 = t\[-den(2) -den(3); 1 0]*t;
	b1 = t\[1; 0];
	c1 = [num(2)-den(2) num(3)-den(3)]*t;
	d1 = 1;
%	[a,b,c,d] = series(a,b,c,d,a1,b1,c1,d1); 
% Next lines perform series connection 
	[ma1,na1] = size(a);
	[ma2,na2] = size(a1);
	a = [a zeros(ma1,na2); b1*c a1];
	b = [b; b1*d];
	c = [d1*c c1];
	d = d1*d;

	i = i + 2;
end

% Take care of any left over unmatched pole pairs.
%   H(s) = 1/(s^2+den(2)s+den(3))
while i < np
	den = real(poly(p(i:i+1)));
	wn = sqrt(prod(abs(p(i:i+1))));
	if wn == 0, wn = 1; end
	t = diag([1 1/wn]);	% Balancing transformation
	a1 = t\[-den(2) -den(3); 1 0]*t;
	b1 = t\[1; 0];
	c1 = [0 1]*t;
	d1 = 0;
%	[a,b,c,d] = series(a,b,c,d,a1,b1,c1,d1);
% Next lines perform series connection 
	[ma1,na1] = size(a);
	[ma2,na2] = size(a1);
	a = [a zeros(ma1,na2); b1*c a1];
	b = [b; b1*d];
	c = [d1*c c1];
	d = d1*d;

	i = i + 2;
end

% Apply gain k:
c = c*k;
d = d*k;

return

function [a,b,c,d] = tf2ss(num, den)
%TF2SS	Transfer function to state-space conversion.
%	[A,B,C,D] = TF2SS(NUM,DEN)  calculates the state-space 
%	representation:
%		.
%		x = Ax + Bu
%		y = Cx + Du
%
%	of the system:
%		        NUM(s) 
%		H(s) = --------
%		        DEN(s)
%
%	from a single input.  Vector DEN must contain the coefficients of
%	the denominator in descending powers of s.  Matrix NUM must 
%	contain the numerator coefficients with as many rows as there are
%	outputs y.  The A,B,C,D matrices are returned in controller 
%	canonical form.  This calculation also works for discrete systems.
%	To avoid confusion when using this function with discrete systems,
%	always use a numerator polynomial that has been padded with zeros
%	to make it the same length as the denominator.  See the User's
%	guide for more details.
%
%	See also: SS2TF.

[mnum,nnum] = size(num);
[mden,n] = size(den);
% Check for null systems
if  (n == 0 & nnum == 0), a=[]; b=[]; c=[]; d=[]; return, end

% Strip leading zeros from denominator
inz = find(den ~= 0);
den = den(inz(1):n);
[mden,n] = size(den);

% Check for proper numerator
if nnum > n
	% Try to strip leading zeros to make proper
	if (all(all(num(:,1:(nnum-n)) == 0)))
		num = num(:,(nnum-n+1):nnum);
		[mnum,nnum] = size(num);
	else
		error('Denominator must be higher or equal order than numerator.');
	end
end

% Pad numerator with leading zeros, to make it have the same number of
% columns as the denominator, and normalize it to den(1)
num = [zeros(mnum,n-nnum) num]./den(1);

% Do the D-matrix first
if length(num)
	d = num(:,1);
else
	d = [];
end

% Handle special constant case:
if n == 1
	a = [];
	b = [];
	c = [];
	return
end

% Now do the rest, starting by normalizing den to den(1),
den = den(2:n) ./ den(1);
a = [-den; eye(n-2,n-1)];
b = eye(n-1,1);
if mnum > 0
	c = num(:,2:n) - num(:,1) * den;
else
	c = [];
end

function [num,den] = zp2tf(z,p,k)
%ZP2TF	Zero-pole to transfer function conversion.
%	[NUM,DEN] = ZP2TF(Z,P,K)  forms the transfer function:
%
%	                NUM(s)
%	        H(s) = -------- 
%	                DEN(s)
%
%	given a set of zero locations in vector Z, a set of pole locations
%	in vector P, and a gain in scalar K.  Vectors NUM and DEN are 
%	returned with numerator and denominator coefficients in descending
%	powers of s.  
%

%	Note: the following will not work if p or z have elements not
%	in complex pairs.

den = real(poly(p(:)));
[md,nd] = size(den);
k = k(:);
[mk,nk] = size(k);
if isempty(z), num = [zeros(mk,nd-1),k]; return; end
[m,n] = size(z);
if mk ~= n
	if m == 1
		error('Z and P must be column vectors.');
	end
	error('K must have as many elements as Z has columns.');
end
for j=1:n
	zj = z(:,j);
	pj = real(poly(zj)*k(j));
	num(j,:) = [zeros(1,nd-length(pj)) pj];
end

function [zd, pd, kd, dd] = bilinear(z, p, k, fs, fp, fp1)
%BILINEAR Bilinear transformation with optional frequency prewarping.
%   [Zd,Pd,Kd] = BILINEAR(Z,P,K,Fs) converts the s-domain transfer
%   function specified by Z, P, and K to a z-transform discrete
%   equivalent obtained from the bilinear transformation:
%
%      H(z) = H(s) |
%                  | s = 2*Fs*(z-1)/(z+1)
%
%   where column vectors Z and P specify the zeros and poles, scalar
%   K specifies the gain, and Fs is the sample frequency in Hz.
%   [NUMd,DENd] = BILINEAR(NUM,DEN,Fs), where NUM and DEN are 
%   row vectors containing numerator and denominator transfer
%   function coefficients, NUM(s)/DEN(s), in descending powers of
%   s, transforms to z-transform coefficients NUMd(z)/DENd(z).
%   [Ad,Bd,Cd,Dd] = BILINEAR(A,B,C,D,Fs) is a state-space version.
%   Each of the above three forms of BILINEAR accepts an optional
%   additional input argument that specifies prewarping. For example,
%   [Zd,Pd,Kd] = BILINEAR(Z,P,K,Fs,Fp) applies prewarping before
%   the bilinear transformation so that the frequency responses
%   before and after mapping match exactly at frequency point Fp
%   (match point Fp is specified in Hz).
%
%   See also IMPINVAR.

%   Gene Franklin, Stanford Univ., motivated the state-space
%   approach to the bilinear transformation.

[mn,nn] = size(z);
[md,nd] = size(p);

if (nd == 1 & nn < 2) & nargout ~= 4	% In zero-pole-gain form
	if mn > md
		error('Numerator cannot be higher order than denominator.')
	end
	if nargin == 5		% Prewarp
		fp = 2*pi*fp;
		fs = fp/tan(fp/fs/2);
	else
		fs = 2*fs;
	end
	z = z(finite(z));	 % Strip infinities from zeros
	pd = (1+p/fs)./(1-p/fs); % Do bilinear transformation
	zd = (1+z/fs)./(1-z/fs);
% real(kd) or just kd?
	kd = (k*prod(fs-z)./prod(fs-p));
	zd = [zd;-ones(length(pd)-length(zd),1)];  % Add extra zeros at -1

elseif (md == 1 & mn == 1) | nargout == 4 %
	if nargout == 4		% State-space case
		a = z; b = p; c = k; d = fs; fs = fp;
		error(abcdchk(a,b,c,d));
		if nargin == 6			% Prewarp
			fp = fp1;		% Decode arguments
			fp = 2*pi*fp;
			fs = fp/tan(fp/fs/2)/2;
		end
	else			% Transfer function case
		if nn > nd
			error('Numerator cannot be higher order than denominator.')
		end
		num = z; den = p;		% Decode arguments
		if nargin == 4			% Prewarp
			fp = fs; fs = k;	% Decode arguments
			fp = 2*pi*fp;
			fs = fp/tan(fp/fs/2)/2;
		else
			fs = k;			% Decode arguments
		end
		% Put num(s)/den(s) in state-space canonical form.  
		[a,b,c,d] = tf2ss(num,den);
	end
	% Now do state-space version of bilinear transformation:
	t = 1/fs;
	r = sqrt(t);
	t1 = eye(size(a)) + a*t/2;
	t2 = eye(size(a)) - a*t/2;
	ad = t2\t1;
	bd = t/r*(t2\b);
	cd = r*c/t2;
	dd = c/t2*b*t/2 + d;
	if nargout == 4
		zd = ad; pd = bd; kd = cd;
	else
		% Convert back to transfer function form:
		p = poly(ad);
		zd = poly(ad-bd*cd)+(dd-1)*p;
		pd = p;
	end
else
	error('First two arguments must have the same orientation.')
end

function [at,bt,ct,dt] = lp2lp(a,b,c,d,wo)
%LP2LP Lowpass to lowpass analog filter transformation.
%   [NUMT,DENT] = LP2LP(NUM,DEN,Wo) transforms the lowpass filter
%   prototype NUM(s)/DEN(s) with unity cutoff frequency to a
%   lowpass filter with cutoff frequency Wo.
%   [AT,BT,CT,DT] = LP2LP(A,B,C,D,Wo) does the same when the
%   filter is described in state-space form.

if nargin == 3		% Transfer function case
        % handle column vector inputs: convert to rows
        if size(a,2) == 1
            a = a(:).';
        end
        if size(b,2) == 1
            b = b(:).';
        end
	% Transform to state-space
	wo = c;
	[a,b,c,d] = tf2ss(a,b);
end

error(abcdchk(a,b,c,d));
[ma,nb] = size(b);
[mc,ma] = size(c);

% Transform lowpass to lowpass
at = wo*a;
bt = wo*b;
ct = c;
dt = d;

if nargin == 3		% Transfer function case
    % Transform back to transfer function
    [z,k] = tzero(at,bt,ct,dt);
    num = k * poly(z);
    den = poly(at);
    at = num;
    bt = den;
end

function [btype,analog,errStr] = iirchk(Wn,varargin)
%IIRCHK  Parameter checking for BUTTER, CHEBY1, CHEBY2, and ELLIP.
%   [btype,analog,errStr] = iirchk(Wn,varargin) returns the 
%   filter type btype (1=lowpass, 2=bandpss, 3=highpass, 4=bandstop)
%   and analog flag analog (0=digital, 1=analog) given the edge
%   frequency Wn (either a one or two element vector) and the
%   optional arguments in varargin.  The variable arguments are 
%   either empty, a one element cell, or a two element cell.
%
%   errStr is empty if no errors are detected; otherwise it contains
%   the error message.  If errStr is not empty, btype and analog
%   are invalid.

errStr = '';

% Define defaults:
analog = 0; % 0=digital, 1=analog
btype = 1;  % 1=lowpass, 2=bandpss, 3=highpass, 4=bandstop

if length(Wn)==1
    btype = 1;  
elseif length(Wn)==2
    btype = 2;
else
    errStr = 'Wn must be a one or two element vector.';
    return
end

if length(varargin)>2
    errStr = 'Too many input arguments.';
    return
end

% Interpret and strip off trailing 's' or 'z' argument:
if length(varargin)>0 
    switch lower(varargin{end})
    case 's'
        analog = 1;
        varargin(end) = [];
    case 'z'
        analog = 0;
        varargin(end) = [];
    otherwise
        if length(varargin) > 1
            errStr = 'Analog flag must be either ''z'' or ''s''.';
            return
        end
    end
end

% At this point, varargin will either be empty, or contain a single
% band type flag.

if length(varargin)==1   % Interpret filter type argument:
    switch lower(varargin{1})
    case 'low'
        btype = 1;
    case 'bandpass'
        btype = 2;
    case 'high'
        btype = 3;
    case 'stop'
        btype = 4;
    otherwise
        if nargin == 2
            errStr = ['Option string must be one of ''high'', ''stop'',' ...
              ' ''low'', ''bandpass'', ''z'' or ''s''.'];
        else  % nargin == 3
            errStr = ['Filter type must be one of ''high'', ''stop'',' ...
              ' ''low'', or ''bandpass''.'];
        end
        return
    end
    switch btype
    case 1
        if length(Wn)~=1
            errStr = 'For the ''low'' filter option, Wn must have 1 element.';
            return
        end
    case 2
        if length(Wn)~=2
            errStr = 'For the ''bandpass'' filter option, Wn must have 2 elements.';
            return
        end
    case 3
        if length(Wn)~=1
            errStr = 'For the ''high'' filter option, Wn must have 1 element.';
            return
        end
    case 4
        if length(Wn)~=2
            errStr = 'For the ''stop'' filter option, Wn must have 2 elements.';
            return
        end
    end
end

